package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class AxisPossionSOTS1D {

	public static void main(String[] args) {
		long begin = System.nanoTime();
		int N = 32;
		double eps = 1.0 / N;
		Mesh1L2 cellMesh = new Mesh1L2();
		cellMesh.line(0.0, 1.0, 128);
		ScalarFunc regionCellFunc1D = (x, llable, param) -> {
			return x[0] < 0.5?1.0:-1.0;
		};
		cellMesh.setDomainLabel(regionCellFunc1D, null, 10);
		FES1L21 fsCell = new FES1L21(cellMesh);
		double[][] FOCS = new double[1][fsCell.GetNdof()];
		double[][] SOCS = new double[1][fsCell.GetNdof()];
		double[][] D = new double[1][1];
		SOTS1DL21 ss = new SOTS1DL21(fsCell, null);
		ScalarFunc Cellcf = (xyz, label, param) -> label == 10? 2: 0.05;
		ss.KxxCellSolution(new ScalarFunc[] {Cellcf}, null, FOCS, D, SOCS);
		
		System.out.println("min: " + MVO.min(FOCS[0]) + "\tmax: " + MVO.max(FOCS[0]));
		System.out.println("Homogenized Kxx: ");
		System.out.println(MVO.toString(D));
		System.out.println("min: " + MVO.min(SOCS[0]) + "\tmax: " + MVO.max(SOCS[0]));
		cellMesh.toTecplot("CellFOCS.dat", FOCS);
		cellMesh.toTecplot("CellSOCS.dat", SOCS);
		
		Mesh1L2 homoMesh = new Mesh1L2();
		homoMesh.line(0, 1, 2048);
		homoMesh.translate(1.0);
		FES1L21 fsHomo = new FES1L21(homoMesh);
		SRMatrix homoA = new SRMatrix(fsHomo.GetNdof());
		fsHomo.assembleStiff(new double[] {D[0][0]}, BVPType.AXISYMMETRIC, homoA);
		double[] homoRHS = new double[fsHomo.GetNdof()];
		fsHomo.assembleSource(new double[] {10.0}, BVPType.AXISYMMETRIC, homoRHS);
		fsHomo.applyBC_MBN(homoA, Direct.X, 2);
		fsHomo.applyBC_MBN(homoRHS, Direct.X, 0.0, 2);
		NewIterSSolver homoSolver = new NewIterSSolver(homoA);
		double[] homoSol = new double[fsHomo.GetNdof()];
		homoSolver.PCGSSOR(homoRHS, homoSol, 1.5, 1);
		
		Mesh1L2 fineMesh1D = new Mesh1L2();
		fineMesh1D.line(0.0, 1.0, 8192);
		fineMesh1D.translate(1.0);
		ScalarFunc regionFunc1D = (x, llable, param) -> {
			double cellX = x[0] * N - (int)(x[0] * N);
			return cellX < 0.5?1.0:-1.0;
		};
		
		fineMesh1D.setDomainLabel(regionFunc1D, null, 100);
		
		ScalarFunc cf = (x, label, param) -> label == 100?2:0.05;
		
		FES1L21 fsFine = new FES1L21(fineMesh1D);
		SRMatrix fineA = new SRMatrix(fsFine.GetNdof());
		fsFine.assembleStiff(new ScalarFunc[] {cf}, null, BVPType.AXISYMMETRIC, fineA);
		
		double[] fineRHS = new double[fsFine.GetNdof()];
		fsFine.assembleSource(new double[] {10.0}, BVPType.AXISYMMETRIC, fineRHS);
		fsFine.applyBC_MBN(fineA, Direct.All, 2);
		fsFine.applyBC_MBN(fineRHS, Direct.X, 0.0, 2);
		
		//fsFine.applyBCZero(fineA, fineRHS, Direct.All, BCType.RCE, 1, 2);
		
		NewIterSSolver solver = new NewIterSSolver(fineA);
		double[] fineSol = new double[fsFine.GetNdof()];
		solver.PCGSSOR(fineRHS, fineSol, 1.5, 1);
		
		fineMesh1D.toTecplot("FineSolution1D.dat", fineSol);
		
		
        //FES2T31 fsFine = new FES2T31(fineMesh);
        double[] HomoFineSol = new double[fsFine.GetNdof()];
        double[] FirstSol = new double[fsFine.GetNdof()];
        double[] SecondSol = new double[fsFine.GetNdof()];
        ss.assembleHeatSOTSSolution(fsFine, HomoFineSol, FirstSol, SecondSol, fsHomo, homoSol, fsCell, FOCS, SOCS, N, eps);
        
        System.out.println("Homo Solution");
        System.out.println("Homo min: " + MVO.min(HomoFineSol) + "\tmax: " + MVO.max(HomoFineSol));
        System.out.println("First Order Approximation: ");
        System.out.println("FOA min: " + MVO.min(FirstSol) + "\tmax: " + MVO.max(FirstSol));
        System.out.println("Second Order Approximation: ");
        System.out.println("SOA min: " + MVO.min(SecondSol) + "\tmax: " + MVO.max(SecondSol));
        
        double[][] wholeSol = new double[4][];
        wholeSol[0] = HomoFineSol;
        wholeSol[1] = FirstSol;
        wholeSol[2] = SecondSol;
        wholeSol[3] = fineSol;
        fineMesh1D.toTecplot("SOTSSolutions.dat", wholeSol);
        double L2Norm = fsFine.computeL2Norm(fineSol);
        double[] error0 = MVO.add(fineSol, -1.0, HomoFineSol);
		double[] error1 = MVO.add(fineSol, -1.0, FirstSol);
		double[] error2 = MVO.add(fineSol, -1.0, SecondSol);
		System.out.println("L2Rrror:");
		System.out.println("Error 0: " + fsFine.computeL2Norm(error0) / L2Norm);
		System.out.println("Error 1: " + fsFine.computeL2Norm(error1) / L2Norm);
		System.out.println("Error 2: " + fsFine.computeL2Norm(error2) / L2Norm);
		
		double H1SemiNorm = fsFine.computeH1SemiNorm(fineSol);
		System.out.println("H1SemiRrror:");
		System.out.println("Error 0: " + fsFine.computeH1SemiNorm(error0) / H1SemiNorm);
		System.out.println("Error 1: " + fsFine.computeH1SemiNorm(error1) / H1SemiNorm);
		System.out.println("Error 2: " + fsFine.computeH1SemiNorm(error2) / H1SemiNorm);		
		long end = System.nanoTime();
		System.out.println("Solve time: " + (end - begin) / 1.0e9 + " s");		
	}

}
